#######################################################
# This script reproduces all R results
# Results are based on R 4.1.3
# Before running this script:
# Set working directory to the project's folder
#######################################################


packages <- c(
  "haven", "lmtest", "arm", "lfe", "foreign", "ggplot2", "viridis",
  "dplyr", "gridExtra", "ggsci", "sp","sf", "RColorBrewer"
)

for (package in packages) {
  if (!require(package, character.only = TRUE)) {
    cat(paste("Package", package, "is not installed.\n"))
    install <- readline(prompt = paste("Do you want to install", package, "? (y/n): "))
    if (tolower(install) == "y") {
      install.packages(package)
      library(package, character.only = TRUE)
    } else {
      cat(paste("Skipping installation of", package, "\n"))
    }
  }
}

source("tools.R")

########################################################
# Map of average Tree Damage by Municipality 
########################################################

vaia <- read_dta("vaia.dta")

load("comunimap")

names(comuni@data)[names(comuni@data) == 'PRO_COM'] <- 'pro_com' 

damage <- merge(comuni, 
                unique(vaia[c("pro_com", "damage_mean", "damage_binary", "damage_max",
                              "damage_fraction", "damaged_prov_binary")]), by = "pro_com")

damaged_prov <- damage[damage@data$damaged_prov_binary == 1 & damage@data$COD_REG != 3, ]

rm(damage)

rm(comuni)

damage_colors <- brewer.pal(9, "Greys")

class_of_damage <- cut(
  damaged_prov@data$damage_mean[damaged_prov@data$damage_binary == 1], 
  
  breaks=quantile(damaged_prov@data$damage_mean[damaged_prov@data$damage_binary == 1],
                  p=c(seq(0,1,length.out=9))), labels=FALSE)

class_of_damage[is.na(class_of_damage)]<-1

pdf(file = "map_damage_mean.pdf")


par(mar=(c(0,0,0,0)))
plot(damaged_prov, col = "white", border = "grey90")

plot(damaged_prov[damaged_prov@data$damage_binary == 1,], 
     
     col = damage_colors[class_of_damage], add = TRUE, border = "grey90")

dev.off()


########################################################
# Coefficient Plots of First Difference Estimates
########################################################

# Incumbent
dt <- read.dta("FD_incumbent.dta")

m1 <- felm(diffshare~ damage_binary|cod_prov|0|SLL_2011, data=dt, subset=anno==2019)

m2 <- felm(diffshare~ damage_binary|cod_prov|0|SLL_2011, data=dt, subset=anno==2014)


m3 <- felm(diffshare~ damage_binary|cod_prov|0|SLL_2011, data=dt, subset=anno==2009)


pdf("FDcoefplot1.pdf", height=5.5/1.9, width=3*6.95/2)

par(mfrow=c(1,3))

par(mar=c(.5,0,.5,0),  tck=-.01)

places <- stackplot(m1, m2,m3, omit=2:24, 
                    suppress.intercept = FALSE, 
                    longnames="" , symmetric = FALSE,
                    force.scale=c(-.03,.06), compressed = TRUE, fifty = FALSE,
                    main="Regional Incumbent: FD",
                    cex.pts = 2, lwd.sd = 4 , cex.main = 1.9
                    
)

text(0.051, places, c("2019","2014", "2009"), cex=1.3 )

# Greens

dt <- read.dta("FD_greens.dta")

m1 <- felm(diffshare~ damage_binary|cod_prov|0|SLL_2011, data=dt, subset=anno==2019)

m2 <- felm(diffshare~ damage_binary|cod_prov|0|SLL_2011, data=dt, subset=anno==2014)

m1g <- m1 

m2g <- m2

m1mat <- cbind( c(coef(m1)),se.coef.felm( m1),c(coef(m2)),se.coef.felm( m2),
                c(NA),c(NA)) 

places <- stackplot(m1mat, omit=2:24, 
                    suppress.intercept = FALSE, 
                    longnames="" , symmetric = FALSE,
                    force.scale=c(-.03,.06), compressed = TRUE, fifty = FALSE,
                    main="Green: FD",
                    cex.pts = 2, lwd.sd = 4 , cex.main = 1.9
                    
)

text(0.02, places[1:2], c("2019","2014"), cex=1.3 )

# M5S
dt <- read.dta("FD_M5S.dta")

m1 <- felm(diffshare~ damage_binary|cod_prov|0|SLL_2011, data=dt, subset=anno==2019)

m1m5s <- m1 

m1mat <- cbind( c(coef(m1)), se.coef.felm( m1),c(NA),c(NA),c(NA),NA ) 

places <- stackplot(m1mat,  
                    suppress.intercept = FALSE, 
                    longnames="" , symmetric = FALSE,
                    force.scale=c(-.03,.06), compressed = TRUE, fifty = FALSE,
                    main="Five Star Movement: FD",
                    cex.pts = 2, lwd.sd = 4 , cex.main = 1.9
                    
)

text(0.02, places[1], c("2019"), cex=1.3 )

dev.off()

########################################################
# Coefficient Plots of Pre-Post DID Estimates (2004-2019 and 2014-2019)
########################################################

d <- read_dta("vaia.dta")

incumbent <- d[d$tipo_elezione == "europea" & 
                 ((d$party == "LN" & (d$cod_reg == 3 | d$cod_reg == 5 | d$cod_reg == 6 | d$cod_prov == 22)) | 
                    (d$party == "SVP" & d$cod_prov == 21)) &
                 d$damaged_prov_binary == 1 & 
                 d$cod_reg != 3,]

m5s <- d[d$tipo_elezione == "europea" & 
           (d$party == "M5S" & (d$cod_reg == 3 | d$cod_reg == 5 | d$cod_reg == 6 | 
                                  d$cod_prov == 22| d$cod_prov == 21)) &
           d$damaged_prov_binary == 1 & d$anno%in%c(2014,2019)&
           d$cod_reg != 3,]

green <- d[d$tipo_elezione == "europea" & 
             (d$party == "green" & (d$cod_reg == 3 | d$cod_reg == 5 | d$cod_reg == 6 | 
                                      d$cod_prov == 22| d$cod_prov == 21)) &
             d$damaged_prov_binary == 1 & d$anno%in%c(2004,2009,2014,2019)&
             d$cod_reg != 3,]

pdf("coefplot1.pdf", height=5.5/1.9, width=3*6.95/2)

dmax.b2 <- felm(party_share ~  st_var(damage_max_t) + as.factor(cod_prov)*as.factor(anno)|pro_com|0|
                  SLL_2011, data=incumbent,cmethod = 'cgm2')

dmean.b2 <- felm(party_share ~  st_var(damage_mean_t) + as.factor(cod_prov)*as.factor(anno)|pro_com|0|
                   SLL_2011, data=incumbent,cmethod = 'cgm2')

dbin.b2 <- felm(party_share ~  st_var(damage_binary_t) + as.factor(cod_prov)*as.factor(anno)|pro_com|0|
                  SLL_2011, data=incumbent,cmethod = 'cgm2')

dmax.b2b <- felm(party_share ~  st_var(damage_max_t) + as.factor(cod_prov)*as.factor(anno)|pro_com|0|
                   SLL_2011, subset=anno>2013, data=incumbent,cmethod = 'cgm2')

dmean.b2b <- felm(party_share ~  st_var(damage_mean_t) + as.factor(cod_prov)*as.factor(anno)|pro_com|0|
                    SLL_2011, subset=anno>2013, data=incumbent,cmethod = 'cgm2')


dbin.b2b <- felm(party_share ~  st_var(damage_binary_t) + as.factor(cod_prov)*as.factor(anno)|pro_com|0|
                   SLL_2011, subset=anno>2013, data=incumbent,cmethod = 'cgm2')

par(mfrow=c(1,3))

par(mar=c(.5,0,.5,0),  tck=-.01)

stackplot( dbin.b2, dbin.b2b, dmax.b2,dmax.b2b, dmean.b2,dmean.b2b, omit=2:24, 
           suppress.intercept = FALSE, 
           longnames="" , symmetric = FALSE,
           force.scale=c(-.03,.03), compressed = TRUE, fifty = FALSE,
           main="Regional Incumbent", pch.type = c(0,1,0,1,0,1),
           cex.pts = 2, lwd.sd = 4 , cex.main = 1.9
           
)

abline(h=c(1.06, 1.14), lty=3)

text(rep(-.028,3), seq(1,1.2,length.out=6),
     c("(Binary)", "Damage", "(Max)",
       "Damage", "(Average)","Damage" ), pos=4, cex=1.5)

dmax.b2 <- felm(party_share ~  st_var(damage_max_t) + as.factor(cod_prov)*as.factor(anno)|pro_com|0|
                  SLL_2011, data=green,cmethod = 'cgm2')

dmean.b2 <- felm(party_share ~  st_var(damage_mean_t) + as.factor(cod_prov)*as.factor(anno)|pro_com|0|
                   SLL_2011, data=green,cmethod = 'cgm2')

dbin.b2 <- felm(party_share ~  st_var(damage_binary_t) + as.factor(cod_prov)*as.factor(anno)|pro_com|0|
                  SLL_2011, data=green,cmethod = 'cgm2')

dmax.b2b <- felm(party_share ~  st_var(damage_max_t) + as.factor(cod_prov)*as.factor(anno)|pro_com|0|
                   SLL_2011, subset=anno>2013, data=green,cmethod = 'cgm2')

dmean.b2b <- felm(party_share ~  st_var(damage_mean_t) + as.factor(cod_prov)*as.factor(anno)|pro_com|0|
                    SLL_2011, subset=anno>2013, data=green,cmethod = 'cgm2')

dbin.b2b <- felm(party_share ~  st_var(damage_binary_t) + as.factor(cod_prov)*as.factor(anno)|pro_com|0|
                   SLL_2011, subset=anno>2013, data=green,cmethod = 'cgm2')

stackplot( dbin.b2, dbin.b2b, dmax.b2, dmax.b2b, dmean.b2,dmean.b2b, omit=2:22, 
           suppress.intercept = FALSE, 
           longnames="" , symmetric = FALSE,
           force.scale=c(-.03,.03), compressed = TRUE,  fifty = FALSE,
           main="Green",
           pch.type = c(0,1,0,1,0,1),
           cex.pts = 2, lwd.sd = 4 , cex.main = 1.9
           
)

abline(h=c(1.06, 1.14), lty=3)

dmax.b2 <- felm(party_share ~  st_var(damage_max_t) + as.factor(cod_prov)*as.factor(anno)|pro_com|0|
                  SLL_2011, data=m5s,cmethod = 'cgm2')

dmean.b2 <- felm(party_share ~  st_var(damage_mean_t) + as.factor(cod_prov)*as.factor(anno)|pro_com|0|
                   SLL_2011, data=m5s,cmethod = 'cgm2')


dbin.b2 <- felm(party_share ~  st_var(damage_binary_t) + as.factor(cod_prov)*as.factor(anno)|pro_com|0|
                  SLL_2011, data=m5s,cmethod = 'cgm2')

stackplot(dbin.b2, dmax.b2, dmean.b2, 
          omit=2:12, 
          suppress.intercept = FALSE, 
          longnames="" , symmetric = FALSE,
          force.scale=c(-.03,.03), 
          compressed = TRUE, main="Five Star Movement",  fifty = FALSE,
          cex.pts = 2, lwd.sd = 4, cex.main = 1.9, pch.type = 1
)
abline(h=c(1.06, 1.14), lty=3)

dev.off()

########################################################
# Event Study Plots for Incumbent and Green
########################################################

d <- read_dta("vaia.dta")

d$damage_binary_04 <- ifelse(d$damage_binary > 0 & d$anno == 2004, 1, 0)
d$damage_binary_09 <- ifelse(d$damage_binary > 0 & d$anno == 2009, 1, 0)
d$damage_binary_19 <- ifelse(d$damage_binary > 0 & d$anno == 2019, 1, 0)
d$damage_mean_04 <- ifelse(d$damage_mean > 0 & d$anno == 2004, d$damage_mean, 0)
d$damage_mean_09 <- ifelse(d$damage_mean > 0 & d$anno == 2009, d$damage_mean, 0)
d$damage_mean_19 <- ifelse(d$damage_mean > 0 & d$anno == 2019, d$damage_mean, 0)
d$damage_max_04 <- ifelse(d$damage_max > 0 & d$anno == 2004, d$damage_max, 0)
d$damage_max_09 <- ifelse(d$damage_max > 0 & d$anno == 2009, d$damage_max, 0)
d$damage_max_19 <- ifelse(d$damage_max > 0 & d$anno == 2019, d$damage_max, 0)

incumbent <- d[d$tipo_elezione == "europea" & 
                 ((d$party == "LN" & (d$cod_reg == 3 | d$cod_reg == 5 | d$cod_reg == 6 | d$cod_prov == 22)) | 
                    (d$party == "SVP" & d$cod_prov == 21)) &
                 d$damaged_prov_binary == 1 & 
                 d$cod_reg != 3,]

green     <- d[d$tipo_elezione == "europea" & 
                 (d$party == "green" & (d$cod_reg == 3 | d$cod_reg == 5 | d$cod_reg == 6 | 
                                          d$cod_prov == 22| d$cod_prov == 21)) &
                 d$damaged_prov_binary == 1 & d$cod_reg != 3,]

results_inc <- combined_results(incumbent) %>%
  mutate(party = "incumbent")

results_green <- combined_results(green) %>%
  mutate(party = "green") %>%
  mutate(measure = ifelse(is.na(estimate), "NA", measure)) %>%
  mutate(estimate = ifelse(is.na(estimate), 0, estimate))

r <- rbind(results_green, results_inc)

r$party <- factor(r$party, levels = c("incumbent", "green"), labels = c("Regional Incumbent", "Green"))

plot <- ggplot(r, aes(x = as.numeric(year), y = estimate, ymin = lower_ci, ymax = upper_ci, shape = measure)) +
  geom_hline(yintercept = 0, color = "black", linewidth=1) +
  geom_vline(xintercept = (2014 + 2019) / 2, linetype = "dashed", color = "gray30", linewidth=1) +
  geom_point(size = 4, position = position_dodge(width = 1.5), stroke=1) +
  geom_linerange(position = position_dodge(width = 1.5), linewidth = 2) +
  labs(x = "", y = "") +
  theme_minimal(base_size = 20) +
  theme(panel.grid.minor = element_blank(),
        panel.grid.major.x = element_blank(),
        panel.grid.major.y = element_line(linetype = "dotted", color="gray30", linewidth = 1),
        legend.position = "bottom",
        strip.text = element_text(size=23, angle = 0, face = "bold")) +
  scale_x_continuous(breaks = c(2004, 2009, 2014, 2019)) +
  facet_wrap(~party, ncol =2) +
  scale_y_continuous(limits = c(-0.007, 0.025)) +
  scale_shape_manual(values = c("damage_binary" = 0, "damage_mean" = 1, "damage_max" = 2, "NA" = 20),
                     labels = c("damage_binary" = "Damage (Binary)", "damage_mean" = "Damage (Average)", "damage_max" = "Damage (Max)", "NA" = "NA")) +
  labs(shape = NULL)

ggsave("combined_ev_study.pdf", plot, width = 15, height = 6, units = "in", dpi = 300)


########################################################
# Event Study Plot of TWFE, IPW-TWFE, and Placebo Estimates
########################################################

d <- read_dta("placebo_vaia.dta")

d$damage_binary <- st_var(d$damage_binary)

d$damage_binary_04 <- ifelse(d$damage_binary > 0 & d$anno == 2004, 1, 0)
d$damage_binary_09 <- ifelse(d$damage_binary > 0 & d$anno == 2009, 1, 0)
d$damage_binary_19 <- ifelse(d$damage_binary > 0 & d$anno == 2019, 1, 0)

# run model without year-province interaction
m.twfe <- felm(party_share ~ damage_binary_04 + damage_binary_09 + damage_binary_19 |pro_com+anno|0|SLL_2011, 
               data = d, cmethod = 'cgm2', weights = NULL)

coefficients <- coef(m.twfe)
conf_intervals <- as.data.frame(confint(m.twfe))

spec  <- rep("twfe", times = 4)
year     <- c("2004", "2009", "2014", "2019")
estimate <- c(coefficients["damage_binary_04"],
              coefficients["damage_binary_09"],
              0,
              coefficients["damage_binary_19"])
lower_ci <- c(conf_intervals[1, 1], conf_intervals[2, 1], NA, conf_intervals[3, 1])
upper_ci <- c(conf_intervals[1, 2], conf_intervals[2, 2], NA, conf_intervals[3, 2])

r.twfe <- data.frame(spec, year, estimate, lower_ci, upper_ci)

# run model with ipweights
d.ipw <- d[!is.na(d$ipw),]

m.ipw <- felm(party_share ~ damage_binary_04 + damage_binary_09 + damage_binary_19 |pro_com+anno|0|SLL_2011,
              data = d.ipw, cmethod = 'cgm2', weights = d.ipw$ipw)

coefficients <- coef(m.ipw)
conf_intervals <- as.data.frame(confint(m.ipw))

spec  <- rep("ipw", times = 4)
year     <- c("2004", "2009", "2014", "2019")
estimate <- c(coefficients["damage_binary_04"],
              coefficients["damage_binary_09"],
              0,
              coefficients["damage_binary_19"])
lower_ci <- c(conf_intervals[1, 1], conf_intervals[2, 1], NA, conf_intervals[3, 1])
upper_ci <- c(conf_intervals[1, 2], conf_intervals[2, 2], NA, conf_intervals[3, 2])

r.ipw <- data.frame(spec, year, estimate, lower_ci, upper_ci)

# run placebo model
d.placebo <- d[!is.na(d$placebo),]
d.placebo$damage_binary_04 <- ifelse(d.placebo$placebo > 0 & d.placebo$anno == 2004, 1, 0)
d.placebo$damage_binary_09 <- ifelse(d.placebo$placebo > 0 & d.placebo$anno == 2009, 1, 0)
d.placebo$damage_binary_19 <- ifelse(d.placebo$placebo > 0 & d.placebo$anno == 2019, 1, 0)

m.placebo <- felm(party_share ~ damage_binary_04 + damage_binary_09 + damage_binary_19 |pro_com+anno|0|SLL_2011,
                  data = d.placebo, cmethod = 'cgm2', weights = NULL)

coefficients <- coef(m.placebo)
conf_intervals <- as.data.frame(confint(m.placebo))

spec  <- rep("placebo", times = 4)
year     <- c("2004", "2009", "2014", "2019")
estimate <- c(coefficients["damage_binary_04"],
              coefficients["damage_binary_09"],
              0,
              coefficients["damage_binary_19"])
lower_ci <- c(conf_intervals[1, 1], conf_intervals[2, 1], NA, conf_intervals[3, 1])
upper_ci <- c(conf_intervals[1, 2], conf_intervals[2, 2], NA, conf_intervals[3, 2])

r.placebo <- data.frame(spec, year, estimate, lower_ci, upper_ci)

# plot
result_combined <- rbind(r.twfe, r.ipw, r.placebo)
result_combined$spec <- factor(result_combined$spec, levels = c("twfe", "ipw", "placebo"))

plot <- ggplot(result_combined, aes(x = as.numeric(year), y = estimate, ymin = lower_ci, ymax = upper_ci, shape = spec)) +
  geom_hline(yintercept = 0, color = "black", linewidth=1) +
  geom_vline(xintercept = (2014 + 2019) / 2, linetype = "dashed", color = "gray30", linewidth=1) +
  geom_point(size = 3, position = position_dodge(width = 1.5), stroke=1) +
  geom_linerange(position = position_dodge(width = 1.5), linewidth=1.5) +
  labs(x = "", y = "") +
  theme_minimal(base_size = 20) +
  theme(panel.grid.minor = element_blank(),
        panel.grid.major.x = element_blank(),
        panel.grid.major.y = element_line(linetype = "dotted", color="gray30", linewidth = 1),
        legend.position = "bottom", 
        plot.title = element_text(face = "bold", size = 20)) +
  scale_x_continuous(breaks = c(2004, 2009, 2014, 2019)) +
  scale_shape_manual(values = c("twfe" = 0, "ipw" = 1, "placebo" = 2),
                     labels = c("twfe" = "TWFE", "ipw" = "IPW", "placebo" = "Placebo")) +
  labs(shape = NULL, title = "Regional Incumbent") 

ggsave("placebo_ev_study.pdf", plot, width = 8, height = 6, units = "in", dpi = 300)

########################################################
# Map of Local Incumbents in Severely Affected Provinces
########################################################

prov19_sf <- st_read("ProvCM01012019/ProvCM01012019_WGS84.shp")

reg19_sf <- st_read("Reg01012019/Reg01012019_WGS84.shp")

# Convert sf objects to SpatialPolygonsDataFrame

prov19 <- as(prov19_sf, "Spatial")

reg19 <- as(reg19_sf, "Spatial")

vaia <- read_dta("vaia.dta")

prov.keep <- vaia[ which(vaia$damaged_prov_binary == 1 & vaia$cod_reg != 3) , ]

prov.keep <- unique(prov.keep[c("cod_prov", "cod_reg")])

aff.prov <- prov19[prov19$COD_PROV%in%prov.keep$cod_prov, ]

aff.prov@data$Incumbent <- ifelse(aff.prov@data$COD_REG == 3, "grey90", 
                                  ifelse(aff.prov@data$COD_PROV == 21, "brown2", "chartreuse3"))

incprov <- factor(aff.prov@data$Incumbent, exclude = NA)

pdf(file = "prov_incumbent.pdf")
plot(reg19[which(reg19$COD_RIP == 1 | reg19$COD_RIP == 2), ], col = "white", border = "grey30")
plot(aff.prov, col= as.character(incprov), border="grey50", add = TRUE)
dev.off()

########################################################
# Plots of Province Trends in Average Vote Share for the Incumbent
########################################################

d <- read_dta("vaia.dta")

incumbent <- d[d$tipo_elezione == "europea" & 
                 ((d$party == "LN" & (d$cod_reg == 3 | d$cod_reg == 5 | d$cod_reg == 6 | d$cod_prov == 22)) | 
                    (d$party == "SVP" & d$cod_prov == 21)) &
                 d$damaged_prov_binary == 1 & 
                 d$cod_reg != 3,]

average_vote_share <- incumbent %>%
  group_by(provincia, anno) %>%
  summarize(average_share = mean(party_share, na.rm = TRUE))

ggplot(average_vote_share, aes(x = anno, y = average_share, group = provincia, color = provincia)) +
  geom_line() +
  geom_point(size = 2) +
  labs(x = "Year", y = "Average Vote Share") +
  scale_x_continuous(breaks = c(2004, 2009, 2014, 2019)) +
  theme_bw() +
  scale_color_lancet() +
  theme(panel.grid.minor = element_blank(),
        axis.ticks.y = element_blank(),
        panel.grid.major.x = element_blank(),
        panel.border = element_blank(), 
        legend.position = "bottom") +
  labs(color = NULL) +
  geom_segment(aes(y = -Inf, yend = -Inf, x = Inf, xend = -Inf), color = "black", linewidth = 0.5)

ggsave(paste0("votebyprovince.png"), width = 8, height = 6, units = "in", dpi = 300)

########################################################
# Plot of Trends in Raw Means of Vote Share by Treatment Group
########################################################

cbPalette <- c("gray60","black")

patterns <- c("incumbent_damage", "green_damage", "m5s_damage")

for (pattern in patterns) {
  
  data_file <- paste0("parallel_", pattern, ".dta")
  png_file <- paste0("parallel_", pattern, ".png")
  
  parallel <- read_dta(data_file)
  trmuni<-factor(parallel$damaged, exclude = NA)
  
  ggplot(parallel, aes(x = year, y = per, colour=trmuni)) +  
    geom_vline(xintercept = (2018+11/12), colour = "gray80", lty = 2, linewidth = 0.8) +
    geom_line(aes(color=trmuni), linewidth = 0.6) + 
    geom_point(size=2) +
    scale_color_manual(name="Municipalities",values=cbPalette,labels=c("Control", "Treated")) +
    scale_x_continuous(name=" ",breaks = c(2004, 2009, 2014, 2019), labels = c("2004", "2009", "2014", "2019"), expand=c(0.035,0.035)) +
    scale_y_continuous(limits = c(0, NA), expand = c(0, 0))+
    ylab("Vote percentage") +
    theme(panel.grid.major.x = element_blank(),
          panel.grid.minor.x = element_blank(),
          panel.grid.major.y = element_line(color = "gray80", linewidth = 0.0025, linetype = 2),
          panel.background = element_blank(),
          axis.line.x = element_line(colour = "black", linewidth = 0.5),
          axis.ticks.x = element_line(colour = "black", linewidth = 0.5),
          axis.ticks.y = element_blank(),
          axis.ticks.length=unit(.2, "cm"),
          axis.title.x = element_text(size = 10), axis.title.y = element_text(size = 10), 
          axis.text.x = element_text(size = 10),axis.text.y = element_text(size = 10),
          legend.title = element_blank(), legend.position = "bottom", legend.key=element_rect(fill="white"),
          plot.title = element_text(face="bold", size=10, hjust = 0.5),
          axis.text=element_text(size=10,color="black"),
          axis.title=element_text(size=12,color="black")) +
    guides(color = guide_legend(override.aes = list(linetype = 0))) + 
    geom_linerange(aes(ymin=cilo,ymax=cihi),linewidth = 0.6) +  
    labs(fill = "")
  
  ggsave(png_file, width = 6, height = 4, units = "in", dpi = 300)
  
  cat("Plot saved for pattern:", pattern, "\n")
}

########################################################
# Event Study Plot Adding Regional Elections
########################################################

d <- read_dta("vaia.dta")

d$damage_binary_04 <- ifelse(d$damage_binary > 0 & (d$anno == 2004 | d$anno == 2005 | d$anno == 2008), 1, 0)
d$damage_binary_09 <- ifelse(d$damage_binary > 0 & (d$anno == 2009 | d$anno == 2010 | d$anno == 2013), 1, 0)
d$damage_binary_19 <- ifelse(d$damage_binary > 0 & (d$anno == 2019 | d$anno == 2020 | d$anno == 2023), 1, 0)
d$damage_mean_04 <- ifelse(d$damage_mean > 0 & (d$anno == 2004 | d$anno == 2005 | d$anno == 2008), d$damage_mean, 0)
d$damage_mean_09 <- ifelse(d$damage_mean > 0 & (d$anno == 2009 | d$anno == 2010 | d$anno == 2013), d$damage_mean, 0)
d$damage_mean_19 <- ifelse(d$damage_mean > 0 & (d$anno == 2019 | d$anno == 2020 | d$anno == 2023), d$damage_mean, 0)
d$damage_max_04 <- ifelse(d$damage_max > 0 & (d$anno == 2004 | d$anno == 2005 | d$anno == 2008), d$damage_max, 0)
d$damage_max_09 <- ifelse(d$damage_max > 0 & (d$anno == 2009 | d$anno == 2010 | d$anno == 2013), d$damage_max, 0)
d$damage_max_19 <- ifelse(d$damage_max > 0 & (d$anno == 2019 | d$anno == 2020 | d$anno == 2023), d$damage_max, 0)

regionali <- d[(d$tipo_elezione == "europea" | d$tipo_elezione == "regionale") & 
                 ((d$party == "LN" & (d$cod_reg == 3 | d$cod_reg == 5 | d$cod_reg == 6 | d$cod_prov == 22)) | 
                    (d$party == "SVP" & d$cod_prov == 21)) &
                 d$damaged_prov_binary == 1 & 
                 d$cod_reg != 3,]

result_combined <- combined_results(regionali)

plot <- ggplot(result_combined, aes(x = as.numeric(year), y = estimate, ymin = lower_ci, ymax = upper_ci, shape = measure)) +
  geom_hline(yintercept = 0, color = "black", linewidth=1) +
  geom_vline(xintercept = (2014 + 2019) / 2, linetype = "dashed", color = "gray30", linewidth=1) +
  geom_point(size = 4, position = position_dodge(width = 1.5), stroke=1) +
  geom_linerange(position = position_dodge(width = 1.5), linewidth = 2) +
  labs(x = "", y = "") +
  theme_minimal(base_size = 20) +
  theme(panel.grid.minor = element_blank(),
        panel.grid.major.x = element_blank(),
        panel.grid.major.y = element_line(linetype = "dotted", color="gray30", linewidth = 1),
        legend.position = "bottom",
        plot.title = element_text(face = "bold", size = 20)) +
  scale_x_continuous(breaks = c(2004, 2009, 2014, 2019), labels = c("Cycle 1", "Cycle 2", "Cycle 3", "Cycle 4")) +
  scale_shape_manual(values = c("damage_binary" = 0, "damage_mean" = 1, "damage_max" = 2),
                     labels = c("damage_binary" = "Damage (Binary)", "damage_mean" = "Damage (Average)", "damage_max" = "Damage (Max)")) +
  labs(shape = NULL, title="Regional Incumbent")

ggsave("regionali_ev_study.pdf", plot, width = 8, height = 6, units = "in", dpi = 300)

########################################################
# Map of Damage Including Lombardia
########################################################

vaia <- read_dta("vaia.dta")

load("comunimap")

names(comuni@data)[names(comuni@data) == 'PRO_COM'] <- 'pro_com' 

damage <- merge(comuni, 
                unique(vaia[c("pro_com", "damage_mean", "damage_binary", "damage_max",
                              "damage_fraction", "damaged_prov_binary")]), by = "pro_com")

damaged_prov <- damage[damage@data$damaged_prov_binary == 1, ]

damage_affected <- damage[damage@data$damage_binary == 1, ]

damage_colors <- rocket(9, d=-1)

class_of_damage <- cut(damage_affected@data$damage_fraction, 9)

damage_colors <- damage_colors[as.numeric(class_of_damage)]

pdf(file = "map_damage_fraction_with_lombardia.pdf", width = 7, height=6)

plot(damage[damage@data$damaged_prov_binary == 1,], col = "white", border = "grey90")
plot(damage_affected, col = damage_colors, add = TRUE, border = "grey90")

dev.off()


